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Recent progress in the development of superconducting circuits has enabled the realization of 
interesting sources of nonclassical radiation at microwave frequencies. Here, we discuss field quadra- 
ture detection schemes for the experimental characterization of itinerant microwave photon fields and 
their entanglement correlations with stationary qubits. In particular, we present joint state tomog- 
raphy methods of a radiation field mode and a two-level system. Including the case of finite quantum 
detection efficiency, we relate measured photon field statistics to generalized quasi-probability dis- 
tributions and statistical moments for one-channel and two-channel detection. We also present 
maximum-likelihood methods to reconstruct density matrices from measured field quadrature his- 
tograms. Our theoretical investigations are supported by the presentation of experimental data, for 
which microwave quantum fields beyond the single-photon and Gaussian level have been prepared 
and reconstructed. 



I. INTRODUCTION 

Microwave frequency quantum fields confined in cavi- 
ties have been generated and characterized with remark- 
able control using Rydberg atoms and superconducting 
qubits for state preparation and readout. These experi- 
ments have illuminated fundamental principles of quan- 
tum physics, e.g. by exploring the coherent superposition 
of quantum states [U [2] and their decoherence [2 O Hj , 
the entanglement between multiple modes [5] and the 
stabilization of Fock states using quantum feedback [B] . 
More recently, progress has been made in the character- 
ization of propagating quantized microwave fields. They 
have so far been prepared in squeezed [7J [3] and single 
photon states [HI [TU] and fully characterized using time- 
correlation measurements [TTJ [T^] and quantum state to- 
mography methods [T3HT7] . These developments have 
also benefited from advances in the efficient detection 
of microwave fields. Both quantum limited linear am- 
plifiers [3 QHH2] and photon counters [23] [24] signifi- 
cantly extend the range of potential quantum optics ex- 
periments using microwave photons interacting with su- 
perconducting qubits, nanomcchanical resonators, quan- 
tum dots [251 ES] , spin ensembles [271 EE] and Rydberg 
atoms [29], 

The use of itinerant microwave photons in quantum 
optics experiments requires efficient field characterization 
methods. A detailed understanding of microwave and op- 
tical field detection schemes allows for adapting existing 
quantum optics tools to the special requirements of mi- 
crowave fields. In the first part of this paper we discuss 
field quadrature detection schemes at microwave frequen- 
cies, their optical analogue and their use for determining 
the quantum state of a single mode of a radiation field. 
We discuss the relation between measurement results of 
single-channel detection schemes and quasi-probability 
distributions. We give new insight into the microwave 
state tomography problem by developing a method to 
reconstruct the maximally-likely Fock space density ma- 
trix directly from the measured probability distributions. 
Furthermore, we present state tomography experiments 
in which quantum states beyond the single photon level 
have been prepared and reconstructed. We also show 



that two-channel microwave detection can be interpreted 
as a positive P function measurement [30] even in the 
presence of added classical detection noise. In the final 
part of the paper we develop new methods, which al- 
low for the characterization of entanglement between a 
localized qubit and a radiation field mode in full joint 
tomography. 

The presented methods are intended for use in state- 
of-the-art circuit QED experiments. However, they are 
also applicable in their general form to other systems 
described by the schematic shown in Fig. [T] It represents 
the generic situation in which two canonically conjugate 
field quadratures X and P of a bosonic mode a and an 
arbitrary spin component a, of a localized qubit are both 
measured. 



II. OPTICAL AND MICROWAVE FREQUENCY 
FIELD DETECTION 

In this section we describe field quadrature detection 
schemes which are frequently used in the optical and in 
the microwave frequency range. We consider both the 
measurement of a single field quadrature and the simul- 
taneous detection of two canonically conjugate quadra- 
tures. We describe the radiation field of interest as a 
single bosonic mode a reaching the detector within a spe- 
cific window in time. The single mode a can be isolated 
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FIG. 1: Schematic drawing of a situation where two conjugate 
field quadatures X and P of a radiation field mode a and an 
arbitrary spin component of a localized qubit a are measured. 
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from the continuum of modes by performing temporal 
mode matching, i.e. integrating the continuous signal 
over the temporal profile of the photon pulse which is 
to be characterized [21 [3T] . We discuss ideal temporal 
mode matching for an exponentially decaying cavity field 
in Appendix [A] 

For a full reconstruction of the quantum state of the 
field both the photon number statistics and all coher- 
ences between the different contributing Fock states have 
to be experimentally determined. This can be achieved 
by measuring generalized field quadrature components 
X$ = \{ae~ l<t> + a^e 1 ^) instead of the photon number 
a* a, which naturally allows for exploring the full phase 
space, i.e. the off-diagonal elements of the density matrix 
in the number state basis |32| . 

In optical systems, where number statistics are natu- 
rally obtained using photon counters, such a field quadra- 
ture measurement can be realized using homodyne detec- 
tion schemes [33] • In this approach, the field of interest 
is combined on a beam-splitter with a strong coherent 
field of a local oscillator, such that the difference of the 
photocurrents at the two beamsplitter outputs are pro- 
portional to a specific field quadrature X^ of the input 
field (see Fig. [2ja)). The quadrature phase <\> can be 
tuned by changing the local oscillator phase. Instead, mi- 
crowave field quadratures are usually measured by down- 




FIG. 2: Field quadrature detection schemes for optical and 
microwave fields, (a) Schematic of balanced optical homodyne 
detection. The signal field a is combined with a coherent local 
oscillator (LO) field with controlled phase <f> at a beamsplit- 
ter and the quadrature amplitude X$ is detected with photon 
counters (n) in the two output arms, (b) Double homodyne 
detection scheme. The signal field a is split into two parts at 
a beamsplitter while introducing an additional vacuum mode 
h. Placing a homodyne detector as described in (a) at each of 
the two beamsplitter outputs allows for measuring two con- 
jugate quadratures (i.e. the complex amplitude S). (c) Mea- 
surement of the complex amplitude at microwave frequencies. 
The signal mode is amplified with a phase-insensitive linear 
amplifier introducing an additional noise mode /lamp- At a 
microwave frequency mixer the amplified output is split into 
two parts, while adding the mode fr m i x , and multiplied with a 
coherent local oscillator field. The down-converted electrical 
field is sampled with analog to digital converters (ADC). 



converting the field with a local oscillator tone using a 
microwave frequency mixer and sampling the electrical 
field directly using analog to digital converters (ADC). 
However, these ADCs are only sensitive enough to de- 
tect large amplitude fields which contain a macroscopic 
number of photons per sampling time, such that a linear 
amplification stage is required in the process of detec- 
tion, as shown in Fig. [2jc) . The noise added during this 
amplification process is typically the main limitation for 
the detection efficiency of microwave fields as discussed 
below. 

Instead of measuring a single field quadrature for dif- 
ferent phases 4>, two conjugate field quadratures can be 
simultaneously measured to get all the information re- 
quired for a complete quantum state reconstruction |3# - 
138] . One possible realization of such a measurement uses 
a beamsplitter and two quadrature detectors at each out- 
put [33] (see Fig.^b)). The beamsplitter necessarily in- 
troduces an additional mode h through its open port. 
This mode adds [at least] the vacuum noise to the signal 
with which the simultaneous detection of conjugate vari- 
ables preserves Heisenberg's uncertainty principle. Tak- 
ing the beamsplitter transformations a — > (a + K)j\[2 
and h (a — h) /y/2 into account, the two detected field 
quadratures at the beamsplitter outputs correspond to 
real X and imaginary P part of the complex amplitude 
a + h) . This holds for both the optical and the microwave 
case. However, for microwaves we still have to consider 
the transformation of the signal mode due to the linear 
amplification stage. A generic phase-insensitive linear 
amplifier transformation can be modeled as jlDHl2^ 

a VGa + y/G - lh{ mp (1) 

where /i amp is an additional bosonic mode accounting for 
the noise added by the amplifier. Again, in the ideal (i.e. 
quantum limited) case /i amp is in the vacuum state, and 
for a more realistic scenario in a thermal state. Combin- 
ing the amplification transformation with the beamsplit- 
ting at the mixing stage (compare Fig.[2|c)) and dividing 
by \J G/2 we find the relation [33] 

S = a + h) = X + iP. (2) 



with the total noise mode h = y ^Q^h amp + J g/lmi x - 
Here, we have defined the complex amplitude opera- 
tor S representing the two conjugate quadratures as 
a single complex number. In the limit of large gain 
G 3> 1 the total noise is dominated by the amplifier 
noise h ss h amp and the following noise contributions can 
be neglected [U]. Furthermore, we notice that once we 
amplify the field phase-insensitively at least the vacuum 
noise is added independently of whether we detect only 
one quadrature or two conjugate quadrature components. 
Once the signal is amplified it is thus natural to detect 
2 conjugate quadratures since the signal-to-noise ratio is 
unaffected by the necessary splitting of the signal. 

It is important to mention that there is a detection 
scheme using linear amplifiers which is ideally noiseless 
for one quadrature component. This is achieved by using 
a phase-sensitive amplifier instead of a phase-insensitive 



3 



one which can, in the quantum limit, be modeled by the 
squeezing transformation [18j HU US] 

a -> VGe^a + (3) 

with the tunable phase <f>. Amplifiers have recently been 
built which are described by this transformation and are 
working close to the quantum limit [7J [55]. The 
quadrature is noiselessly amplified while its conju- 
gate quadrature is deamplified. The detection scheme 
is thus equivalent to an optical homodyne detection 

pi on [mug. 

We note that while for optical fields the simultaneous 
detection of two conjugate quadratures requires a more 
complicated setup than for photon number detection it is 
the natural measurement observable for microwave fields 
which we will therefore focus on in this work in the con- 
text of quantum state reconstruction. 

III. QUANTUM STATE RECONSTRUCTION 
BASED ON SINGLE CHANNEL FIELD 
QUADRATURE DETECTION 

Here, we describe quantum state tomography based on 
the measurement of the complex field amplitude S. The 
goal of quantum state reconstruction is the estimation 
of the density matrix p a which characterizes the state 
of the field mode a. This is experimentally achieved by 
preparing the state many times and performing a set of 
measurements on these states, which contain information 
about the diagonal and off-diagonal elements of p a ■ De- 
pending on the set of observables the obtained results 
have a direct analogy to particular representations of the 
density matrix. In the case of field quadrature detection 
these representations are phase space distributions such 
as the Husimi-Q function or the Wigner function [37H39"] . 
in the following we discuss their relation to statistical mo- 
ments and the Fock basis representation of the density 
matrix. 



A. Phase space distributions 

Due to the non-orthogonality of coherent states 
(a|/3) = e - sH 12 1 1 e Q ^ an arbitrary density matrix 
p a can be expanded as a linear combination of projectors 
onto coherent states 

Pa= f P a (a)\a)(a\. (4) 

J a 

Here we have defined / = J c d 2 a for integrals over the 
complex plane and P a (a) as the Glauber-Sudarshan P 
function [35J HO] > which uniquely represents the density 
matrix as a distribution in phase space. P a (ct) is always 
real- valued but can be negative and can contain singular- 
ities proportional to derivatives of the Dirac S distribu- 
tion to all orders [35]. As can be seen from its definition 
Eq. Q the P function reduces to a two-dimensional Dirac 
distribution P a (a) — S^(a — /3) for coherent states |/3). 
Coherent states thus appear as single points in phase 



space with no statistical spread similar to their classi- 
cal counterparts. For this reason and due to its possible 
negative values the P function does not directly describe 
the statistics of measurements. However, it is very useful 
since its statistical moments directly correspond to the 
normally ordered moments 

((a^) m a n ) = f (a*) m a n P a (a) (5) 

J a 

of the field operator a and because of its analogy to prob- 
ability distributions of classical fields. A second distribu- 
tion, which is of particular relevance for the following 
discussion, is the Husimi-Q function 

Q a (a) = -(a\p a \a), (6) 

since it generates the anti-normally ordered moments 

(aVD = f (a*) m a n Q a (a) . (7) 

J a 

Substituting Eq. Q into the definition of the Q function, 
we note that it is related to the P function by a Gaus- 
sian convolution. For coherent states Q a (a) becomes a 
two-dimensional Gaussian distribution with variance 1 
centered around the coherent state amplitude. Half of 
these fluctuations describe the intrinsic vacuum fluctu- 
ations of the quantum field, the other half describe the 
minimal added uncertainty when directly measuring a Q 
function, which requires the simultaneous detection of 
two conjugate field quadratures. 

Both distributions are special cases of the s- 
parametrized quasi-probability distribution (QPD) 
W a (a,s) 

Q a (a)=W a (a,-l) (8) 
P a (a) = W a (a, 1). (9) 

which has been introduced by Cahill and Glauber [37] as 
a generalized phase space representation of the density 
matrix where the parameter s 6 (— oo,+l]. For differ- 
ent values of s the QPDs are related to each other by a 
Gaussian convolution [37] 

W a (a,s) = g^^cxp (~ 2| "~f |2 ) W a W,t) (10) 
where t > s. 

An intuitive interpretation of the parameter s relates 
to the amount of fluctuations which are contained in the 
distribution in units of half photons. For s = 0we obtain 
the Wigner function W a (a) = W a (a, 0) which include 
the intrinsic vacuum fluctuations but no additional noise 
due to measurement. In the case s = 1 we identify the 
P function where even the vacuum fluctuations are not 
represented. On the other hand, for s = — 1 the QPD 
corresponds to the Q function where both the vacuum 
fluctuations and the minimal added detection noise are 
embedded. As discussed below, additional classical de- 
tection noise leads to s < — 1 when identifying measured 
distributions with a generalized QPD. 
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B. Measurement of phase space distributions 

In order to understand the relation between general- 
ized QPDs and measured distributions in different exper- 
imental situations we assume that the complex amplitude 
S = a + iv as introduced above is the measured observ- 
able and that the results S of repeated measurements 
are stored in a 2-dimensional histogram D^(S) where 
the two axes represent the real and imaginary part of S. 
From this measured distributions all statistical moments 
can be numerically evaluated as 



{S*) n S m D^(S) . (11) 



If the noise added by the detection chain is independent 
of the signal generated by the photon source the signal 
mode a and the noise mode h are uncorrelated p = p a <£> 
Ph- Under this assumption the moments of the measured 
distribution can be decomposed into products of signal 
and noise moments 

rn,n / \ / \ 

i(&) n s m ) Pa = (JJ (jJ<(^W)<a m ~V) n ~ j ), 

(12) 

Here, we have chosen an operator ordering where the 
signal moments (a m (a ji ) n ) appear anti-normally ordered 
and the noise moments ((h)) m h n ) normally ordered. 
Note that since S is a normal operator [S, £ft] — one 
can express Eq. ( 12 ) also with opposite ordering as shown 
in Eq. (16) later in the text. 



The probability distribution for the sum of two inde- 
pendent random variables a + h) is identical to the con- 
volution of the individual distributions for a and h*. As 
a result, one possible representation of the probability 
distribution (S) is given by the convolution [51] 



D W (S)= / Ph(S* -a*)Q a (a). 



(13) 



In the following we discuss special cases of Eq. (13). At 



optical frequencies the measurement of S can be realized 
using a double homodyne or heterodyne detection and 
the noise mode h is nearly in the vacuum state for which 
P h {fi) = 8W{fi) resulting in 



(14) 



In contrast, for microwave fields the noise mode h is often 
in a thermal state with mean photon number Nq ranging 
typically from 0.5 to 10 if parametric or SQUID amplifiers 
are used [H 123 E2] or between 30 an 200 if the first 
amplification is performed by a transistor based amplifier 

Q31H3H1HS]. In this Case Ph ^ = e ~ lo,l2/N °/nN acts 
as a Gaussian filter and by comparing with Eq. ( 10 ) we 
obtain the broadened QPD 



D^\S) = W a (S,-l-2N Q ). 



(15) 



Note that finite thermal noise in h can be equivalently 
interpreted as optical homodyne detection with finite de- 
tection efficiency r\ for which the measured distribution 



is given by D^(S) = W a (S, l-lr]' 1 ) ^M- Added noise 
can thus be understood as a reduced detection efficiency 
?7 = l/(l + ivo). 

We conclude that under the experimentally verified as- 
sumption [141 115) of h being in a thermal state not cor- 
related with a the measured distribution of S is a direct 
measurement of the generalized QPD and therefore con- 
tains all information required to reconstruct the density 
matrix p a of the state of interest or to test its nonclassi- 
cal properties [MjEjjj. In contrast to other reconstruction 
schemes only a single observable S needs to be measured. 

However, in many experiments the mean photon num- 
ber of the noise field is larger than the mean photon 
number of the signal field No > (a' a) and consequently 
the features of measured probability distributions are on 
first sight dominated by the noise distribution. There- 
fore the goal is to systematically extract the information 
contained about mode a in the measured QPD and rep- 
resent it in a form, which allows for a direct estimation 
of the properties of the state, such as the fidelity with 
respect to an expected density matrix. 



C. Determination of normally ordered moments 



One way of quantifying the properties of a quantum 



of 



state is to analyze the statistical moments ((a')" 
the field operator [T3 EE]> since quantities such as the 
mean amplitude, the mean photon number and the vari- 
ance in the photon number can be extracted immedi- 
ately. In this section we discuss the approach developed 
in Ref. |14j to extract these moments from the measured 
distributions in the presence of significant amplifier noise 
N . The basic idea is to deconvolve the QPDs for the 
field operators a and h order by order. 

Rewriting Eq. (12) with a different choice of operator 
ordering 



((&rs m ) Pa = mm<(« t M<fc n ~W TO_i >> 

(16) 

we find that once the anti-normally ordered moments of 
the noise mode (h n (W) m ) are known, the set of linear 
equations can be solved for ((a') n a m ). From Eq. ( 11 ) we 



note that a reference measurement D^°^°^(S), for which 
a is prepared in the vacuum, gives direct access to the 
moments {h n Qv) m ), since all normally ordered moments 
in a with n,m ^ are then ((a^) n a m ) — and Eq. (16) 
reduces to 



>|0>(0| 



(17) 



In cryogenic setups such a reference measurement with a 
in the vacuum can typically be performed by cooling the 
source of radiation into the ground state or very close to 
it [ST]. The identity in Eq. (17) can be understood as 
follows: The situation with a in the vacuum state corre- 
sponds to an ideal Q function measurement for the noise 
mode h and the moments generated by this distribution 
are exactly the anti-normally ordered ones appearing in 
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Eq. ( 17 ) . We finally invert Eq. ( 16 ) to extract the desired 
moments ((a') n a m ) of the mode to be characterized. 

In principle, the moments of the measured histograms 
can be evaluated to arbitrary order. However, there 
are limitations in the accuracy with which the moments 



((a 



can be determined depending on the integra- 



tion time and the detection efficiency As investigated 
theoretically in Ref. the statistical error of the mo- 
ments increases with increasing order. The result shows 
that the number of measurements which are necessary 
to extract a moment of order M with a given precision 
scales with (1 + JVo) . The measurement time necessary 
to determine higher order moments with a fixed precision 
thus scales exponentially with increasing order. 

The state of a single mode of the radiation field has 
an infinite number of degrees of freedom, i.e. an infinite 
dimensional Hilbert space. This makes it in principle im- 
possible to exactly reconstruct a state, because an infinite 
amount of information is to be acquired. However, the 
measurement of a finite set of moments often allows for 
a controlled reduction of the relevant state space [56]. 



D. Special classes of states and the Fock space 
density matrix 

One class of states which is characterized by a finite set 
of moments comprises coherent, thermal and squeezed 
(i.e. Gaussian) states, for which the statistical moments 
up to second order 



{{a), (J a), (a*)} 



(18) 



determine all higher order moments. In order to analyze 
how close the reconstructed state really is to a Gaussian, 
one has to measure the third order cumulants and eval- 
uate their deviations from zero. 

A second class of states which can be reconstructed us- 
ing a finite set of measured moments includes those with 
finite photon number contributions satisfying (n\p a \m) = 
for m,n > N in the Fock basis {|n)}. For these states 
the normally ordered moments 



((a f ) n a m ) = to or n>N 



vanish and the state is completely determined by the fi- 
nite set of moments 



{((a t )"a m )} to and n<N. 



(20) 



It is important to note that it necessarily follows from 
((cr) N a N ) — that there are no Fock states \n) with n > 
N contributing to the density matrix. If ((a') N a N ) < e 
can be verified experimentally one knows an upper bound 



n>N ^ '' n>N 

for the sum of higher order Fock state populations. The 
approximation made when truncating the Hilbert space 
is thus well-controlled. 



If such a truncation is possible the moments can be 
mapped to a density matrix in Fock representation by 
evaluating the sum [58] 



(m\p a \n) = 



1 



£ 



\J n\m\ /! 

M{{(a)) n a m )) 



(22) 



up to terms of order 2A^. 

The described procedure is very efficient since the eval- 
uation of moments from the measured distributions as 



well as finding the solution of Eq. ( 16 ) requires only small 



(19) 152] 



computational effort. Furthermore the moment represen- 
tation provides a very intuitive picture to extract funda- 
mental properties of the quantum state. 



IV. MAXIMUM-LIKELIHOOD STATE 
ESTIMATION BASED ON GENERALIZED 
COMPLEX AMPLITUDE DETECTION 

Due to the unavoidable statistical imprecision in ex- 
pectation values extracted from a finite number of mea- 
surements, a direct mapping from the measurement data 
to the desired state representation does not in general re- 
sult in a completely positive density matrix. Maximum- 
likelihood state estimation aims to correct for that. In 
this section we discuss two different maximum-likelihood 
procedures applicable to complex amplitude detection 
schemes as relevant for the circuit QED experiments un- 
der consideration. The first method is based on the ex- 
perimentally determined finite set of moments {(cv) n a m ) 
together with their respective standard deviations <5 n , m . 
The second one estimates the density matrix directly 
from the measured probability distributions. 



Maximum-likelihood procedure based on 
measured moments 



In order to find the most likely density matrix given a 
set of measured moments and their respective standard 
deviations b n ^ m we maximize the log-likelihood function 



£ Log = -£ Ta-K^rO - Tr[p a (at)"a™]| 2 (23) 



S 2 



with respect to the elements of the density matrix p a . 
The properties p a > and Trp a — 1 of the density ma- 
trix are included as constraints in the maximization of 
Eq. (23). The standard deviations 5 n , m appear in the 
denominator of each term, such that moments which 
are determined with low accuracy contribute to the log- 
likelihood function with less weight. 

This maximization problem can be formulated as a 
semi-definite program, for which efficient numerical so- 
lutions exist jini ISO] • Note that this maximum likelihood 
scheme is particularly efficient for states which contain 
only few photons since in this case only a finite set of 
moments is non-zero. 
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We have tested the described maximum-likelihood pro- 
cedure based on experimental data sets obtained in a 
circuit QED experiment. In addition to the generation 
of single photon states [U] we have prepared two pho- 
ton Fock states and their coherent superposition with 
the vacuum. Note that for the reconstruction of 2 pho- 
ton states it is necessary to measure photon correlations 
including moments up to sixth order. The accurate mea- 
surement of ((a^) 3 a 3 ) - compared to previous measure- 
ments in which ((a^) 2 a 2 ) had been the highest order mea- 
sured moment [TTJ [T^l [2] - was enabled by a Josephson 
parametric amplifier used as the first amplifier in the de- 
tection chain |HT] . 

Based on the measured moments and their respective 
standard deviations up to order n + m = 8 we reconstruct 



each density matrix by maximizing Eq. (23). In order to 



demonstrate that higher order photon number popula- 
tions are not relevant for the description of the state if 
one of the diagonal moments (a)) N a N ) is measured to 
be close to zero (compare Eq.(21|), we have chosen a 
Hilbert space with up to four photon Fock states. The 
results (see Fig. [3k) show that only the zero, one and 



two photon Fock states contribute to the reconstructed 
density matrices while the higher Fock states stay un- 
populated. A compromise between the size of the Hilbert 
space and the likelihood of the reconstructed state may 
be found by applying the Akaike or Bayesian information 
criterion [55] to reduce the complexity of the model used 
for reconstructing the state. In order to illustrate the 
quantum character of the reconstructed states we have 
transformed the density matrices into their correspond- 
ing Wigner functions [14] , which show negative values in 
all four cases (see Fig.[3j3). 

We estimate the statistical error in the fidelities of 
the reconstructed density matrices by repeating the 
likelihood maximization for resampled sets of moments 
[STfl, [S3] . The resulting standard deviations of the resam- 
pled fidelities are below 2% for all reconstructed states. 
The small statistical errors are due to the high overall 
microwave detection efficiency of r\ = 0.19 of our setup 
in combination with the large number of measurements 
exceeding 10 s for all the shown density matrices. Since 
high repetition rates of up to 10 MHz [11] are possible 
for circuit QED experiments we believe that the maxi- 
mum likelihood approach is well suited in this context. 
However, in experiments for which only a small number 
of samples is available alternative methods such as the 
Bayesian approach |64l 165] may be advantageous com- 
pared to maximum likelihood procedures. 



B. Iterative maximum-likelihood procedure based 
on measured histograms 

In addition to the moments based maximum likelihood 
scheme we formulate an iterative procedure which esti- 
mates the density matrix directly from the measured his- 
tograms. This reconstruction method is useful for photon 
states which contain a large number of contributing Fock 
states and consequently a large number of non- vanishing 
moments. In addition to this practical relevance it gives 



insight into the interpretation of the measured probabil- 
ity distribution. 

The measurement of a quantum observable can be de- 
scribed by a set of positive operator valued measures 
(POVM) lij [SS], which have the property that the prob- 
ability pj for getting the respective measurement result 
is given by pj = Tr[pTlj]. The operators Tlj need to be 
positive and hermitian but not necessarily projectors. In 
the ideal case they form a decomposition of the Hilbert 
space JV Hj — t. Preparing and measuring a system in 
state p repeatedly will return each of the possible results 
fj times. The most likely state pml given this set of data 
is the one which maximizes the likelihood function 



£ = nTr[pn i ] 



(24) 



Note that in order to find a unique global maximum of 
C, it is a necessary condition that an arbitrary density 
matrix can be constructed as a linear combination of H, . 
As a counterexample, if the POVM are given by a com- 
plete orthogonal set of projectors fij = the ML 
function C is independent of the off-diagonal elements of 
p expressed in the \j) basis. The maximization of C can 
thus only identify the most likely diagonal density matrix 
elements (j\p\j)- 

It is computationally demanding to directly determine 
Pml for high-dimensional Hilbert spaces. However, the 
density matrix /?ml carL be found using iterative meth- 
ods [351 [57-69 . In order to formulate the ML iteration 
procedure we define the operator 



Rip) = E 



fj 



Tr[pH 



(25) 



The iterative method for updating the density matrix 
[S3I7D] 



Pk+i = AfR(p k )pkR(Pk) 



(26) 



with renormalization constant A/", has shown good con- 
vergence towards /?ml [22] • As an initial condition for 
the iteration procedure one either chooses the maximally 
mixed state po — where d is the dimension of the 
reconstructed Hilbert space, or constructs a more realis- 
tic initial condition by taking into account the measured 
moments. 

In a practical implementation where the phase space 
is discretized and the Hilbert space is truncated to finite 
dimensions we might also be faced with the situation that 
the POVM operators do not sum to the identity operator 
TV Ilj = G ^ 1 . In this situation the iteration procedure 
can be modified as 



Pk+1 =NG- 1 R{p k ) Pk R{ Pk )G- 



(27) 



to guarantee convergence towards the most likely density 
matrix [Til 172"] . 



1. Iterative method for ideal complex amplitude detection 

The method described above has been adapted to op- 
tical homodync detection by Lvovsky [68 in 2004 and is 
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FIG. 3: (a) Absolute value of the experimentally reconstructed density matrices (grayscale) in comparison with the ideal 
ones (wireframes) for the four indicated quantum states, (b) Measured density matrices transformed into their corresponding 
Wigner functions W(x,p). The fidelities between ideal states |i/>) and measured density matrices are evaluated as F = {ip\p\il>). 



frequently used in experiments based on optical homo- 
dyne tomography [75H75] , Here we adapt the method to 
measurements of the complex amplitude operator S. We 
start with the case of ideal detection, i.e. for the noise 
mode h being in the vacuum state. 



As discussed in Section III B the measured probability 
distribution in this case is the Q function D^(S) = 
Q a (S). The underlying set of POVMs n s is thus defined 
by the condition 



Q a (S) =Tr[p a n s ] 



(28) 



Since the Q function can be written as the expectation 
value Q a {a) = ^(a\p a \a) with respect to coherent states 
| a) we identify the well-known result [75] 



n 



S=a 



|a)(a| 



(29) 



Here and in the following we have labeled the possible 
measurement results of S by a to emphasize their relation 
to coherent states. 

The coherent state projectors fl Q have both the de- 
sired properties: They sum up to the identity operator 
J H a = 1 and they allow for the construction of an arbi- 
trary density matrix as a linear combination of projectors 



Pa 



P a (a)U c 



(30) 



compare with Eq. Q. Based on this knowledge we can 
directly apply the iteration procedure Eq. (27). 



Full state tomography thus requires the measurement 
of only a single observable S which ideally projects onto 
coherent states. Due to the properties of coherent states 
all information about the phase of the field necessary to 
reconstruct the off-diagonal density matrix elements is 



contained in this measurement. This is one of the reasons 
why the discussed detection scheme has great potential 
in microwave photon field tomography - especially since 
the advent of nearly quantum-limited amplifiers [7J IT5I 



2. Iterative method for generalized complex amplitude 
detection 



Due to noise added by amplifiers as well as finite ra- 
diation losses in waveguides and microwave components 
the mode h is typically not described by the vacuum but 
a thermal state with mean photon number N$. In the 
following we show how to reconstruct the density matrix 
p a in this situation. We keep the discussion as general as 
possible and allow for mode h being in an arbitrary state 
described by ph which can be specified experimentally 
using a reference measurement. The following procedure 
has to the best of our knowledge not yet been discussed 
in literature. 

Preparing the signal mode a in the vacuum state we 
can measure the Q function of mode h since Z)[|o)(o|] ( a ) _ 
Qhioi*). Applying the iterative maximum likelihood 
scheme for ideal detection we reconstruct the most likely 
state for the noise mode ph- To account for this noise 
state in the reconstruction of p a we identify the modified 
POVM operators ITa h ', which describe the measurement 
process under the condition that the detection system is 
in state ph- The result, which can be shown by verifying 
the identity 



Tr[ Pa n[? hl ] =£>W(a) 



Eq.Jl3) 



(31) 



between POVMs and the expected measured distribu- 
tion, is 



n[f 1 = -T h (a)p h Tl(a). 



(32) 



Here we have defined the displacement operator Th(a) = 
e ah ~ a h and ph as the most likely density matrix with 
respect to the reflected histogram Qh{— a*). Note that 
since displaced Fock states are orthonormal and complete 



T{a)\m)(n\T\a) = U n , 



(33) 



the relation J a Ida = 1 holds for any valid detector 
state ph- This leads to the remarkable result that the 
reconstruction method is unbiased for arbitrary detector 
states. The only two requirements for the method to ap- 
ply are that the signal and noise modes are uncorrelated 
P = Pa ® Ph and that o can be cooled into the vacuum 
state or any other known state. Both of these conditions 
can be realized experimentally to good approximation as 
discussed before. In order to estimate the density matrix 
p a we can again apply the iterative method using fla h ^ 
as a set of POVMs. 

We have applied the iterative maximum likelihood 
scheme to the same data sets as presented in Fig(3] and 
found quantitative agreement between the two methods 
to about 1%. As described in the following, we have 
tested the iterative method also for a coherent state |a) 
with mean amplitude a rs 1.7, for which we expect higher 
photon number states to be occupied. We first apply 
the iterative procedure to the reference histogram which 
characterizes the detector state ph- Its diagonal elements 
Pn = ( n \Ph\n) are shown in Fig.|4]as dots which are very 
well described by a thermal distribution (solid line) with 
mean photon number No ~ 4.4. The off-diagonal ele- 
ments (not shown) are all smaller than e = 0.004. There- 
fore, the detection noise is very well approximated by 




0.05 



FIG. 4: Diagonal matrix elements p n = {n\ph\n) of the de- 
tector state (dots) obtained with the iterative maximum like- 
lihood method from experimental data. The photon number 
distribution is well described by a thermal distribution (solid 
line) with mean photon number iVo ~ 4.4. The inset shows 
the reconstructed density matrix with fidelity F = 95% of a 
coherent state with a ~ 1.7. 



thermal noise. Taking into account the estimated detec- 
tor state ph we construct Ida and iterate the maximum 
likelihood procedure for the coherent state histogram. 
The resulting estimated density matrix p a is shown in 
the inset of Fig. [4] and has a fidelity of F = 95% com- 
pared to an ideal coherent state. 

Note that in order to reconstruct and express the den- 
sity matrix of the detector state ph with high accuracy 
we have to take into account a Hilbert space of a dimen- 
sion which is approximately 10 times the noise number 
Nq. It is therefore numerically challenging to implement 
the iterative procedure in cases where the noise number 
is large. If this is the case one should preferably work 
with the moments based maximum likelihood method 
presented in the previous section. 



TWO CHANNEL DETECTION AND THE 
POSITIVE P DISTRIBUTION 



We have already pointed out that field quadrature 
measurement is the most commonly used detection 
method for microwave frequency fields. Due to its well- 
established implementation it has also been possible to 
experimentally realize Hanbury Brown Twiss-type setups 
(see Fig. [5]) where two instead of one complex ampli- 
tudes are measured [ill H2 H3 EE] [79] . The advantage of 
such a detection scheme is that ideally the system noise 
in the two detection arms is uncorrelated and only the 
signal mode a contributes to the cross-correlations be- 
tween the two output arms. In this section we provide 
a quantum optics description of a generic two-channel 
microwave detection chain [T5] [73] as shown in Fig. [5] 
We formulate the main advantages of such a measure- 
ment setup and show that under reasonable assumptions 
a direct measurement of a positive P distribution |80j is 
realized. This relation between the positive P distribu- 
tion and the two-channel detection scheme gives impor- 
tant insight into the general statistical properties of the 
obtained measurement results. 




FIG. 5: Two channel detector with radiation incident from 
input mode a. Each of the beamsplitter outputs has an indi- 
vidual amplification stage adding noise in modes hi and fe. 
In both channels the complex amplitude is measured [111143] . 
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A. Two-channel detection 

The main difference between the one- and the two- 
channel setup depicted in Fig. [5] is the additional beam- 
splitter which splits the signal mode a equally into two 
parts while introducing an additional mode v. As a re- 
sult, the input modes at the two amplifiers are given by 
(a ± v)/y/2 and the total measured complex amplitudes 
Si and S2 can be expressed as 



51 = a + v + y/2h\, 

5 2 = a — v + V2h,2, 



(34) 



where hi and h 2 are the modes accounting for the sys- 
tem noise of the two detection chains. Under the as- 
sumptions justified below that (i) the mode v is in the 
vacuum state, that (ii) all other modes are uncorrelated 
P = Pa® Phi ® phi-, an d that (iii) the noise has no phase- 
coherence (h™) = = (h™), Vm > we find the follow- 
ing cross-correlations [35] 



((SjrS*™) = ((a t ) m a n ). 



(35) 



The above assumptions require that (i) the open beam- 
splitter port is connected to a bath of zero temperature, 
that (ii) the signal is not correlated with the two com- 
pletely independent amplifier chains, and that (iii) the 
noise does not depend on the phase defined by the refer- 
ence local oscillator, all of which can in good approxima- 
tion be experimentally realized [TTJ [T2] 

This result is remarkable since under realistic condi- 
tions the above cross-correlations completely describe the 
state of mode a without any influence of the detector 
noise modes. This means that the scheme is largely in- 
dependent of the choice of amplifiers and noise sources, 
even if the noise constitutes the majority of the power 
in the signals Si , S2 ■ As we show in the following these 
properties can be understood in terms of the positive P 
function representation [80j . 



B. The positive P function 

The positive P function was introduced by Drummond 
and Gardiner [50] as a theoretical concept for the solution 
of Fokker-Planck equations. In contrast to the Wigner 
function and Glauber-Sudarshan P function it is com- 
pletely positive and has all properties of a genuine prob- 
ability distribution. The positive P function P(a,/3*) is 
defined as a non-diagonal expansion of the density matrix 
in the coherent state basis 



Pa = / P(a,f3)\ 



(36) 



Like the P function it generates the normally ordered 
moments of the field operator 



((at)"V>=/ P(a,/3)a"(/ry 

Ja,0 



(37) 



Furthermore this four-dimensional probability distribu- 
tion P(a,f3*) can be shown to be positive, not unique 
and to exist for any quantum state |80) . To resolve the 
problem of uniqueness one can resort to the canonical 
choice |35j of the positive P function which is given by 



47T 



" B -" 9 W^V (38) 



While the positive P function is often considered artifi- 
cial and only of theoretical relevance, Braunstein et ai. 
have shown that it can be interpreted as the probability 
distribution for the simultaneous measurement of four 
quadrature variables |35j . A scheme for an optical exper- 
iment was proposed by Agarwal [3D] based on fourfold 
balanced homodyne detection. To our knowledge this 
scheme has so far not been implemented, probably due 
to the significant experimental effort necessary. 



C. Two channel detection as a measurement of the 
positive P function 

The scheme by Agarwal and the two channel mi- 
crowave detection scheme in Fig. [5] are equivalent up 
to the presence of the amplifier noise. Furthermore un- 
der the assumptions made above about the noise, the 
observables Si,S 2 generate the normally ordered mo- 
ments in the same way as the positive P function. It 
is thus natural to assume that the probability distribu- 
tion of the measurement data P(S\, S2) is a positive P- 
representation of the input mode a. In appendix [B] we 
calculate this distribution and find 



P(Si, S2 



Si* - p* 



s 2 * - p* 
~7T 



(39) 



where Qi, 2 {ct) are the Q functions of the system noise 
modes hi,h 2 . When the noise added to both channels is 



in a thermal state with mean photon number Nq Eq. ( 39 ) 
simplifies to 



P(Si,S 2 ) 



exp ( 

in (N 



\S!-S2\' 

4(jV +l) 



W a 



Si + s 2 



2N 
(40) 



and for quantum limited detection, i.e. Nq = 0, to 



P(Si,S 2 ) 



1 

-in 



exp 



I Si — S2 



Qa 



Si + s 2 



-Fcan (Si, S2) 



(41) 
(42) 



The compelling result is that for quantum limited de- 
tection the measurement data distribution corresponds 
to the canonical choice of the positive P-representation 



P can (a, /?). Moreover, we show in Appendix B 2 that the 



measured distribution is always a positive P distribution 



1 



P(Si,S 2 ) 



Si,S 2 



HSi)(S 2 |a) 
(S 2 |Si) 



?a(«) (43) 
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for any thermal populations Ni , N2 in the detector noise 
modes. As a consequence, the density matrix can be 
directly evaluated from the measured P(S\,S2) using 
Eq. ( 36 ) even in the presence of significant thermal noise 



of unequal powers in the detection chains. These results 
suggest that the measurement of a positive P distribution 
is possible with current microwave frequency quadrature 
detection setups such as the one described in Ref . [TT] . 



VI. JOINT TOMOGRAPHY SCHEME FOR A 
QUBIT - PHOTON FIELD SYSTEM 

In the previous sections we have described experimen- 
tal schemes for characterizing microwave radiation on 
a quantum level. In many experiments where such ra- 
diation fields are relevant they interact with stationary 
quantum objects such as superconducting qubits. Simi- 
lar to optical photons in atomic systems [8TH87] . itiner- 
ant quantum microwave fields have the potential to act 
as a quantum channel to connect spatially separated su- 
perconducting circuits with each other. In this section 
we discuss a state tomography scheme to reconstruct the 
joint density matrix of the qubit-photon field system [55] 
using linear detection. The scheme is an extension of the 
photon field tomography methods described in the pre- 
vious sections and is applicable to any system in which 
the complex field amplitude S and the qubit state pop- 
ulation along an arbitrary axis can both be measured in 
each trial of an experimental run. The presented method 
allows to include finite detection efficiencies for both the 
qubit and the photon field detection. 



A. Qubit state tomography 

In order to describe the joint tomography scheme we 
first discuss the concepts of qubit tomography. For the re- 
construction of the qubit density matrix p a one measures 
the Pauli expectation values (at) along the different spin 
axes Oi £ {<J x ,<j y ,a z } which in the measurement basis 
{|0 Z ),|1 Z }} are represented by the corresponding Pauli 
matrices. After state preparation the qubit is rotated 
such that the desired spin component ai points along the 
measurement axis. This rotation is followed by a read- 
out procedure during which the measurement result is 
encoded in a classical quantity q [89 . In the context of 
circuit QED single shot read-out [55l [5U] is not always 
available and q can take a continuous spectrum of val- 
ues where depending on the qubit state each value has a 
probability Di{q) to occur. Here the index i £ {x, y, z} 
specifies the measurement basis. 

The distribution Di(q) obtained after repeating the 
measurement many times can be fitted to the weighted 
sum of 2 reference distributions po (q) for the ground and 
Pi (q) for the excited state 

A (?) = l ~ 1 Y A Pi (?) + 1± ^Po (?) (44) 

to extract the Pauli expectation values (<Tj). Based on 
these values the density matrix can then be determined 



as p a = |(1 + X)i( cr i) cr i)- Note that instead of using 
Eq. ( 44 1 the qubit population can also be extracted from 
the mean values of q [91] as done in many experiments . 



B. Joint tomography 

A joint tomography scheme is expected to allow for a 
full characterization of the system also when photon field 
and qubit are correlated with each other. The goal is to 
determine all matrix elements (s, n\p a>a \s', m) of the joint 
density matrix p a ^ a where s, s' £ {0 Z , 1 2 } label the qubit 
basis states and n, m the photon field number states. It 
can be shown that these matrix elements are uniquely 
determined by the set of moments ((cv) n a m ai) in the 
following way 



(0z,m\p <7ya \0 z ,n} 
(l z ,m\p a;a \l z ,ri) 
(lz,m|/9 o ., o |0 z ,n) 



-M(((^ra m ) + ((a^ra m a z )) 
\M (((^"a" 1 ) - ((at)«a'V z » 
~M(((aira m (a x +ia v )}) (45) 



Here, M. is the linear map from the moments to the 
density matrix in the number state basis as defined in 
Eq. ( 22 1 . The scheme described in the following allows 
to measure all the necessary moments in Eq. ( 45 1 . 



We consider the case that in each trial of an exper- 
imental run both q, characterizing the qubit state, and 
the complex amplitude S, characterizing the photon field, 
are measured. For each state preparation both numbers 
are stored in a 3D histogram Di(S,q) where the index 
i labels the chosen qubit rotation before measurement. 
To evaluate the desired expectation values we first deter- 
mine the Pauli expectation values (<Jj)s conditioned on 
the complex amplitude result S. This is done by fitting 
each trace of Di (S, q) along the q axes to the calibration 
histograms po(q) and p\{q). Based on the knowledge of 
(cj)s we determine the photon field distributions condi- 
tioned on a specific qubit measurement result as 



D U (S) = M LJpk. D i (S, q) (46) 



where Ao and A/i are appropriate normalization con- 
stants which guarantee that J S D .(S) = J s D li (S) = 1. 
For example, Dq^ (S) is the photon field distribution un- 
der the condition that the qubit is measured with result 
in the x-basis. 

Given these histograms one can evaluate the condi- 
tioned moments ({&) n S m )\0i and ((S r ) n S m )\U using 
Eq. (11). If signal field and qubit are not correlated 
with the noise p a ^ a ® ph we can use the techniques de- 
scribed in section |III| to extract the desired quantities 
(((r) n a m <Ji). In recent experiments we have used the 
presented method to analyze correlations between single 
itinerant microwave photons entangled with a supercon- 
ducting qubit [BTj . 
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C. Maximum-likelihood state estimation for a joint 
system 

As for photon state tomography, it is desirable to find 
a set of POVM operators for the described measurement 
scheme. This allows to construct an iterative maximum- 
likelihood state estimation procedure and furthermore 
provides insight into the conditioned histograms intro- 
duced in Eq. ((46 1. 



For a perfect single shot qubit readout with a binary 
measurement result or 1, the POVMs are given by 
the projectors onto eigenstates of the Pauli operators, 
fto; = |0j)(0j| and = |lj)(lj|. They are complete in 
the sense that an arbitrary qubit density matrix can be 
explicitly written as a linear combination of projectors 

/v = |(i + £>i)(n 0i -nij). 

Including the photon field measurement the total set 
of POVMs is 



n. 



(47) 



with s € {0, 1} and i £ {cc, y, z} for which the expectation 
values with respect to the total density matrix are related 
to the measured histograms by 



Tr[p CT , a n Qi0 J 
Tr[p (7 „.n„ 



1 



J cr,a lx a ,li 



(48) 



with (<7j) being the unconditioned Pauli expectation val- 
ues. Note that this remains valid for qubit readout with 
limited single shot fidelity since the storage of the data 
in 3-dimensional histograms allows for capturing all the 
necessary qubit-photon correlations. By fitting the 3D 
histogram data along the g-axes to the expected ground 
and excited state distributions (see Eq. ( 44 ) ) we account 



for the finite read-out efficiency. Using the set of POVMs 
given in Eq. ( 47 1 we are able to use the iterative max- 



imum likelihood procedure described above to estimate 
the most likely density matrix for the combined system. 



VII. CONCLUSIONS 
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Appendix A: Temporal mode matching 

Here we discuss the relation between the single time- 
independent mode a describing the photon pulse which 
is to be characterized by state tomography and the time- 
dependent field clout (i) which is continuously sampled in 
an experiment. The link between the two is given by a 
mode-matching relation 



a = / dtf(t)a ou t(t) 



(Al) 



where the normalization condition J d£|/(i)| 2 = 1 of the 
temporal profile function f(t) guarantees that [a, a*] = 1 
is satisfied. The best choice of f(t) depends on the tem- 
poral shape of the field, i.e. the properties of the coupling 
between radiation source and the bath under observa- 
tion. In the following we discuss optimal temporal mode 
matching for a single-sided cavity acting as the radiation 
source. 

We represent the cavity mode by the creation operator 
A(t), which in a Heisenberg picture is time-dependent. 
We assume that at time t = the cavity is prepared in a 
specific state described by the statistics of 4(0) and then 
left under free evolution 92 

4(i) = e~^4(0) + ^e _i r / dre^a in (r). (A2) 



From input-output theory [93] we know that the cavity 
field decays with rate k into the output modes according 
to 



In summary, we have presented schemes for the char- 
acterization of a single radiation field mode and its en- 
tanglement with a two-level system based on linear am- 
plification and quadrature detection. We have discussed 
single channel and two-channel detection schemes and 
showed that the latter enables a direct measurement of 
a positive P function even in the presence of significant 
added amplifier noise. For both the photon field and the 
joint tomography scheme we have discussed maximum- 
likelihood procedures which take into account the full 
measured quasi-probability distributions. 

Due to the recent progress in the development of 
quantum-limited amplifiers and microwave photon coun- 
ters we believe that the investigation of itinerant mi- 
crowave photons has a great potential in quantum sci- 
ence. The state tomography methods described in this 
paper enable the experimental characterization of such 
microwave radiation fields and of optical fields detected 
with finite efficiency. 



CW (i) = VftA(*) - <*«(*)• ( A3 ) 

The input modes a- m (t) can be understood as a continu- 
ous stream of independent modes each reaching the res- 
onator at time t and ideally carrying only the vacuum 
noise. 

By inserting the above expressions into the definition 
of a we get only one term 4(0) yj~k J °° dte~ Kt depending 
on the cavity field. In order to maximize the efficiency 
in detecting the state prepared at time t = we have to 
find f{t) which maximizes this term. The choice fit) = 
\f~HeJ~ ^9(t) does so, where O(i) is the Heaviside step 
function. The total expression then reduces to 

/■OO /"OO 

a = 4(0)re / dte-^-n 1 / 2 e-^a in it)dt 

ft Jo (A4) 

+ k 3/2 / e~ Kt / e^tHnOrJdrd*. 
Jo Jo 
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which due to the identity 



poo pi 

K 3/2 / e -«i / e T ain ( T )drdi 
Jo Jo 

poo / poo \ 
= KV2 J \J @{ J- T ^ Ktdt ) e¥a in(T)dT (A5) 



k 1 / 2 / e '» a;nfr)dr. 



simplifies to 



a = A(0). 



(A6) 



By proper choice of /(t) we can thus recover the state 
of the source field A(0) with unit efficiency in the trans- 
mission line. Note that a finite mode matching efficiency 
only reduces the total detection efficiency but does not 
affect the statistical properties of a. 



Substituting this equation and the integral forms of the 
characteristic functions in the definition of P(S\, S 2 ) we 
get 

P(S U S 2 )= [ P«G9)Qi(»?)Q2(7)A (B8) 

where 

D = 7r~ 4 f exp(2i(v^ry + 0* - S£) - c.c.)) 



exp( Z2 (^2 7 + /3*-5*)-c.c.)). 



(B9) 



which reduces Eq. (|B8| to 



P(Sx,S 2 ) 



1 



Pa(/?)Ql 



s* - p* 



Q-2 



s 2 * - 



(BIO) 



Appendix B: Probability distribution for two 
channel complex envelopes 

Here we calculate the joint probability distribution of 
the complex envelopes in a two channel detection scheme 
along the lines of Ref. (3U] ■ By definition the probability 
distribution of the measurement data Si , S 2 is given by 
the Fourier transform of the characteristic function 



P(Si,S 2 ) 



1 



Z 1 Sl+Z 2 S2—Z\S 1 —Z 2 S 2 



where 



Xss(zi,z 2 ) 



_ / ziS\+z 2 §l -z{Si-z 2 S2 



Xss(zi,z 2 ). 

(Bl) 



(B2) 



By substituting the operator Eqs. (34 1 for Si and S 2 we 
find 

Xss(zi,z 2 ) = (B3) 
Xa(zi + z 2 )\v(zi - z 2 )xh 1 (-%/2zi)xh 2 (-%/2z2) 

where we introduced the characteristic functions for the 
four different modes as 



Xa(z) 

Xv(z) 

Xh t { Z ) 



P a {(i)e 



z/3*-z*p 



z* hi 



me 



z/3*-z*/3 



(B4) 

(B5) 
(B6) 



We can simplify these expressions by introducing the 
following physical assumptions. First, mode v is assumed 
to be in the vacuum state. In this case its characteristic 
function is the identity and we have 



Xss{zi,z 2 ) = Xa(zi +z 2 )x/ ! . 1 (-V2zi)x/ l2 (-V / 2^2)- 



(B7) 



1. Thermal noise 

In the next step we assume that the noise modes h%, h 2 
are in thermal states e'^ /Ni+1 /n{Ni + 1) 94 with 
mean photon numbers Ni , N 2 . The probability distri- 
bution of the measurement data is then 



P(S 1 ,S 2 )= / P a {p) 



exp 



( \Si-p\- 

\ 2(JVi + l 



\S 2 -0\ 2 
) 2(N 2 + 1) 



47r 2 (7Vi + l)(A^ 2 + l) 

(Bll) 

Comparing this with the formula for the s-parametrized 
QPD 

W a (S,s) = ^ / F Q (/3)exp f- 2 -\^l) (B12) 



1 - s 



we can identify the relation 



1 - s 



P(S 1 ,S 2 ) = -^e- lJ i^W a (S,s) (B13) 



2TTN tl 



by defining 



Ni + l N 2 + 1 Q 



s= -1- 



4 7Vi7V 2 + 2iVi + 2iV 2 



Ntot =N!+N 2 + 2. 



(B14) 

(B15) 
(B16) 



If we have the same noise level on both channels N\ 
N 2 = N Q , we find s = -1 - 2N and thus 

S1-S2I 2 
e 4 < N o+i) /Si + So 

(BIT) 

In the case of quantum limited detection Ni = N 2 = 
we have s = —1 and the distribution of our mea- 
surement data corresponds to the canonical positive P- 
representation of mode a 

Si + S 2 



I isi-s 2 i 2 
P(S 1 ,S 2 ) = —e- lJL ^Q a 



(B18) 
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2. Equivalence to canonical positive P function 



To proof that Eq. (B13l is also a positive P function 
when the thermal noise levels are unequal N\ ^ N 2 we 
show that 35 



- f P(S 1 ,S 2 ) 
77 J Si,S 3 



(g|g 1 )(5 2 |g) 
(Sal Si) 



Q{a). (B19) 



Using Eq. (Bill, the lhs of Eq. (B19) is 



CX P (-fe+S - SkSl) (a\Si)(S 2 \a) 



US u s 2 PaW) ^(N 1 + l)(N 2 + l) (S 2 \Si 



This is a multi-dimensional Gaussian integral in the vari- 
ables Si , 5*2 and can be solved to give 



-|/3-a| 2 



(B20) which is exactly the Q function. 



(B21) 
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